Revise implementation of ht_neuron #491
Conversation
…s code; some code tidying.
# Conflicts: # models/ht_neuron.cpp
# Conflicts: # models/ht_connection.h # models/ht_neuron.cpp # models/ht_neuron.h
…example parameter value.
| , tau_P_( 50.0 ) | ||
| , delta_P_( 0.2 ) | ||
| , tau_P_( 500.0 ) | ||
| , delta_P_( 0.5 ) |
There was a problem hiding this comment.
Since this could potentially have an influence on existing user scripts, could you justify the change in default parameters?
There was a problem hiding this comment.
The actual values of these two parameters are not given in the Hill-Tononi (2005) paper. The new values are the default values in the original Synthesis simulator code by Sean Hill. I sent a mail to NEST User on 22 September warning of changes to the Hill-Tononi model and asking for comments from anyone who may be affected, but have not received a single answer, so I think we are safe here.
| @@ -49,8 +49,7 @@ RecordablesMap< ht_neuron >::create() | |||
| Name( "Theta" ), &ht_neuron::get_y_elem_< ht_neuron::State_::THETA > ); | |||
| insert_( | |||
| Name( "g_AMPA" ), &ht_neuron::get_y_elem_< ht_neuron::State_::G_AMPA > ); | |||
| insert_( | |||
| Name( "g_NMDA" ), &ht_neuron::get_y_elem_< ht_neuron::State_::G_NMDA > ); | |||
| insert_( Name( "g_NMDA" ), &ht_neuron::get_g_NMDA_ ); | |||
There was a problem hiding this comment.
Why is the recordable g_NMDA different from the state variable? (cf. get_g_NMDA)
EDIT: ok, I understand, the state variable G_NMDA is the g(t) from the docstring and the recordable is the real g_NMDA = g(t)*m(V, t), which cannot be stored directly as a state variable since it is not what's updated.
But then the naming looks misleading to me, I guess we should change G_NMDA to something else and update the docstring... though I must say I'm at a loss regarding the name. Maybe G_T would be the most sensible choice...
There was a problem hiding this comment.
Naming is a bit of a mess here. In the paper, it is \tilde{g}_NMDA = m(V) g_NMDA(t), where g_NMDA(t) is the double exponential function.
I think that the end user will be most interested in the actual conductance of NMDA channels when Mg-blocking is taken into account (the m(V)). Therefore, I think the recordable should be g_NMDA and the getter get_g_NMDA_(), but I think I will change the name of the state variable, maybe to G_NMDA_TIMECOURSE.
There was a problem hiding this comment.
Yes, I agree that g_NMDA should be the variable recorded, and G_NMDA_TIMECOURSE looks good!
| * NMDA conductance | ||
| * | ||
| * We need to take care to handle instantaneous blocking correctly. | ||
| * If the unblock-variables Mg_{fast,slow} > Mg_ss, the steady-state |
There was a problem hiding this comment.
This looks unclear to me, you mean "if the unblock variables are greater than their respective steady-state values", right? I think the combination of equation and text sentence with commas is rather disturbing.
There was a problem hiding this comment.
You are right, fixed.
| , NMDA_Vact( -58.0 ) // mV | ||
| , NMDA_Sact( 2.5 ) // mV | ||
| , NMDA_Tau_1( 4.0 ) // ms | ||
| , NMDA_Tau_2( 40.0 ) // ms |
There was a problem hiding this comment.
If you think it's worth it (I'm not sure), having all taus in lowercase would be more coherent, though it would break existing codes... and Tau_m instead of tau_m is really too bad since it's supposed to be common to many models...
| @@ -280,7 +300,10 @@ class ht_neuron : public Archiving_Node | |||
| DG_GABA_A, | |||
There was a problem hiding this comment.
I would be in favour of changing G_NMDA to G_T but there should at least be a comment like
G_NMDA, // state variable corresponding to g(t), with g_NMDA = m(V, t) * g(t)
or something equivalent.
| @@ -409,6 +432,26 @@ class ht_neuron : public Archiving_Node | |||
| { | |||
| return S_.I_h_; | |||
| } | |||
| double_t | |||
There was a problem hiding this comment.
I suppose because the branch is a bit old .... Fixed.
|
@Silmathoron I now noticed that the model needs a bit more work. I will change all names to use But I probably won't get that done today. |
|
Ok, thanks for your hard work! |
| * We need to take care to handle instantaneous blocking correctly. | ||
| * If the unblock-variables Mg_{fast,slow} > Mg_ss, the steady-state | ||
| * value for the present membrane potential, we cannot change those values | ||
| * in State_[], since the ODE Solver may call this function multiple times |
There was a problem hiding this comment.
I would prefer it if "those values" were specified
There was a problem hiding this comment.
I will fix this with the next update to the PR.
|
The documentation of the NMDA current is very clear and helpful! Although I would add an explicit definition of g (what is g_NMDA in the paper). The calculation of the NMDA part of I_syn and propagation of Mg_fast and Mg_slow look correct to me! I approve this pull request. |
|
|
||
| m(V) = 1 / ( 1 + exp( - ( V - NMDA_Vact ) / NMDA_Sact ) ) | ||
| where g(t) is a beta function (difference of exponentials). | ||
|
|
There was a problem hiding this comment.
I would prefer it if g was stated explicitly
There was a problem hiding this comment.
I will fix this with the next update to the PR.
# Conflicts: # nestkernel/nest_names.cpp # nestkernel/nest_names.h
# Conflicts: # nestkernel/nest_names.cpp # nestkernel/nest_names.h
|
@Silmathoron I have now completed my revisions, so you can continue reviewing. |
There was a problem hiding this comment.
Nice restructuring!
Here is a first round of comments, I'll try to check the equations and test the model next week:
- in
ht_connection.h, line 44, it should betau_P(ortau_rec, if possible, cf. comments below). - the
THIS MODEL NEURON HAS NOT BEEN TESTED EXTENSIVELY!comment leaves me perplex... what do you mean exactly? Wouldn't a test file solve the problem?
| def< double >( d, "tau_P", tau_P_ ); | ||
| def< double >( d, "delta_P", delta_P_ ); | ||
| def< double >( d, "P", p_ ); | ||
| def< double >( d, names::tau_P, tau_P_ ); |
There was a problem hiding this comment.
Since we are trying to limit the number of names and given that this is clearly a recovery time, I think it would be a good idea to set it to names::tau_rec and specify in the docstring "the recovery time tau_P in the paper is called tau_rec here", or something equivalent
There was a problem hiding this comment.
I see your point, but since everything here is about P, and I would like to stay reasonably close to the paper to ease comparisons, I would like to stick with tau_P.
| def< double >( d, names::theta_eq, theta_eq ); | ||
| def< double >( d, names::tau_theta, tau_theta ); | ||
| def< double >( d, names::tau_spike, tau_spike ); | ||
| def< double >( d, names::spike_duration, t_spike ); |
There was a problem hiding this comment.
t_spike is not detailed anywhere, could you add a line of explanation in the docstring?
| //! Timer (counter) for potassium current. | ||
| int r_potassium_; | ||
| //! Timer (counter) for spike-activated repolarizing potassium current. | ||
| int r_spike_; | ||
|
|
||
| bool g_spike_; //!< active / not active |
There was a problem hiding this comment.
The name g_spike_ seems weird for a boolean telling whether the K current is in action or not. Also this variable is not documented though it's a recordable...
| - Intrinsic currents I_h (pacemaker), I_T (low-threshold calcium), | ||
| I_Na(p) (persistent sodium), and I_KNa (depolarization-activated | ||
| potassium). | ||
| - Several apparent typographical errors in the descriptions of |
There was a problem hiding this comment.
I guess you mean errors in the paper [1], but maybe it should be more explicit...
| t_peak = tau_2 tau_1 / ( tau_2 - tau_1 ) ln( tau_2 / tau_1 ) | ||
|
|
||
| m(V, t) = a(V) m_fast*(V, t) + ( 1 - a(V) ) m_slow*(V, t) | ||
| a(V) = 0.51 - 0.0028 V |
There was a problem hiding this comment.
There is an apparent problem with the homogeneity, here... I suppose 0.0028 is in mV^{-1}, but could you specify it explicitely? (maybe introduce a named constant).
Shouldn't the user want to be able to set these two values? If not (which would be surprising to me), would you mind explaining their origin?
There was a problem hiding this comment.
I sure should document this better and provide details in a notebook. Given the number of parameters in this model that are fixed as numerical constants, making them all settable would make the code more complex (just inventing good names for all constants would be a major effort). Therefore, I would like to leave them as they are (seems to be quite usual for these more complex models). The particular values here are from Vargas-Caballero and Robinson, J Neurophysiol 89: 2778–2783, 2003. I will add the reference.
| @@ -155,6 +175,8 @@ const Name growth_rate( "growth_rate" ); | |||
| const Name gsl_error_tol( "gsl_error_tol" ); | |||
|
|
|||
| const Name h( "h" ); | |||
| const Name h_E_rev( "h_E_rev" ); | |||
| @@ -179,6 +207,9 @@ const Name Interpol_Order( "Interpol_Order" ); | |||
| const Name interval( "interval" ); | |||
| const Name is_refractory( "is_refractory" ); | |||
|
|
|||
| const Name KNa_E_rev( "KNa_E_rev" ); | |||
| @@ -202,8 +233,18 @@ const Name N_channels( "N_channels" ); | |||
| const Name n_events( "n_events" ); | |||
| const Name n_proc( "n_proc" ); | |||
| const Name n_receptors( "n_receptors" ); | |||
| const Name NaP_E_rev( "NaP_E_rev" ); | |||
| const Name needs_prelim_update( "needs_prelim_update" ); | ||
| const Name neuron( "neuron" ); | ||
| const Name NMDA_E_rev( "NMDA_E_rev" ); |
| @@ -288,6 +332,8 @@ const Name synapse_model( "synapse_model" ); | |||
| const Name synapse_modelid( "synapse_modelid" ); | |||
| const Name synaptic_elements( "synaptic_elements" ); | |||
|
|
|||
| const Name T_E_rev( "T_E_rev" ); | |||
|
All should be fixed now. I am continuing work on the notebook, next step is a systematic test of the intrinsic currents (voltage-clamping is coming to NEST ...) |
I_H verified in notebook.
…better clarity and proper initialization with steady-state values.
… channels tested successfully and decent similarity to Fig 5 in paper achieved.
…HT_NMDA # Conflicts: # nestkernel/nest_names.cpp # nestkernel/nest_names.h
|
@Silmathoron @ingablundell I have just pushed a major revision of the model including a notebook exploring all aspects of the model (docs/model_details). Enjoy! Note that failing tests may be due to Travis trouble this afternoon. |
There was a problem hiding this comment.
Great work on that, the update function is much clearer now!
There are only a few points left that need some clarifications; also I realized that I forgot to ask, but I find it very surprising that, in such a detailed model, conductances cannot be taken directly from biological data... could you at least provide typical values that would allow for an approximate conversion?
| "\\end{equation}\n", | ||
| "This equation is also given in a comment in Synthesis, but is missing from the implementation.\n", | ||
| "\n", | ||
| "**Note: NEST implements the equation according to Compte et al (2003) with $N_{NaP}=3$, while Synthesis uses $N_{NaP}=1**\n", |
There was a problem hiding this comment.
Missing $ closure after $N_{NaP}=1
| "\n", | ||
| "#### Synaptic \"minis\"\n", | ||
| "\n", | ||
| "Synaptic \"minis\" due to spontaneous release of neurotransmitter quanta [HT05, p 1679] are not included in the NEST implementation of the Hill-Tononi model. This due to the fact that the total mini input rate for a cell was just 2 Hz and they cause PSP changes by $0.5 \\pm 0.25$mV only and thus should have minimal effect.\n", |
There was a problem hiding this comment.
For people working on cultures (like me) this is a very important information, so maybe it would be good to make it appear clearly at the beginning of the notebook.
(also minis are usually proportional to the neuron's in-degree, so this constant value is a bit surprising to me... could you specify the conditions?)
There was a problem hiding this comment.
From the paper: "In addition, we modeled “minis” ... as low-amplitude PSPs (mean 0.5+-0.25 mV) .... The mean frequency of these Poisson distributed synaptic minis was set to 2 Hz (total for an individual cell)."
| * function multiple times and in arbitrary temporal order. We thus need to | ||
| * use local variables for the values at the current time, and check the | ||
| * state variables once the ODE solver has completed the time step. | ||
| * If the unblock-variables m_NMDA_{fast,slow} are greater than the |
There was a problem hiding this comment.
Naming is not in agreement with lines 87--88 (m_{fast,slow}_NMDA)
|
|
||
| const double D_influx = | ||
| D_influx_peak / ( 1.0 + std::exp( -( V - D_thresh ) / D_slope ) ); | ||
| return P_.tau_D_KNa * D_influx + D_eq_0; |
There was a problem hiding this comment.
D_eq_0 was not present before and the value returned by the function is apparently not homogeneous... could you explain the role of this new variable?
There was a problem hiding this comment.
Renamed D_eq as in the paper/equations.
There was a problem hiding this comment.
What about the homogeneity?
I noticed that, indeed, the equation for D does not seem to satisfy this requirement, unless D_influx has a dimension of
There was a problem hiding this comment.
I hadn't thought about this, but you are right. The dimension of D does not matter, as long as it is the same as of d_1/2, since the two are divided in the end, so dimensions drop out. I suggest to leave D and d_1/2 dimensionless, since they represent some "factor" which has no immediate physical counterpart (but models some sort of ion concentration). Then D_eq also is dimensionless, while D_influx and D_influx,peak must have units of inverse time, explicitly 1/ms. This makes sense, since these values essentially represent ion concentration change. I have updated the notebook.
| def< double >( d, names::theta, y_[ THETA ] ); // Threshold | ||
| } | ||
|
|
||
| void | ||
| nest::ht_neuron::State_::set( const DictionaryDatum& d, const Parameters_& ) | ||
| nest::ht_neuron::State_::set( const DictionaryDatum& d, | ||
| const ht_neuron& node, |
There was a problem hiding this comment.
I don't understand this modification of the function
There was a problem hiding this comment.
Several of the intrinsic channels have slow dynamics, O(1s) for I_h, e.g. If you change parameters, especially the membrane potential using SetStatus/SetDefaults, it will often be useful to set the state variables for these channels to their new steady-state values. Otherwise, one would have to simulate several seconds of biological time to get to equilibrium. Especially when setting model parameters at the beginning of a simulation, the new equilibration capability should be useful. When changing parameters during the simulation, on the other hand, you probably want to see how the system responds to the change (eg my voltage stepping tests and the wake-sleep transitions in HT05), so then one will not want to equilibrate.
There was a problem hiding this comment.
I'm sorry, my previous statement was unclear! I don't understand why you added const ht_neuron& node as parameter, could you not simply call the functions?
There was a problem hiding this comment.
set() is a member function of State_, not of ht_neuron. Therefore, I cannot directly call member functions of ht_neuron. But since I now call methods in ht_neuron to get the equilibrium values, I don't need the Parameters_ explicitly and will remove that argument.
There was a problem hiding this comment.
Of course, I read too fast, sorry...
| @@ -577,7 +624,7 @@ nest::ht_neuron::Buffers_::Buffers_( const Buffers_&, ht_neuron& n ) | |||
| nest::ht_neuron::ht_neuron() | |||
| : Archiving_Node() | |||
| , P_() | |||
| , S_( P_ ) | |||
| , S_( *this, P_ ) | |||
There was a problem hiding this comment.
... and same answer, essentially: the state should be initialized to the correct steady-state values.
There was a problem hiding this comment.
Here I need both parameters, since I compute some equilibrium values explicitly in the constructor.
| {E_rev,g_peak}_{h,T,NaP,KNa} - reversal potential and peak conductance for | ||
| intrinsic currents | ||
| receptor_types - dictionary mapping synapse names to ports on | ||
| neuron model | ||
| recordables - list of recordable quantities. | ||
| recordables - list of recordable quantities | ||
| equilibrate - if given and true, time-dependent activation |
There was a problem hiding this comment.
Maybe a precision that values are otherwise initialized at zero would be good
There was a problem hiding this comment.
I added and "otherwise" part.
| @@ -285,12 +288,15 @@ class ht_neuron : public Archiving_Node | |||
|
|
|||
| double g_peak_KNa; | |||
| double E_rev_KNa; // mV | |||
| double tau_D_KNa; // ms; in P_ for technical reasons, not meant to be public | |||
There was a problem hiding this comment.
I turned it into a normal, accessible parameter now.
|
I have just pushed a new version addressing @Silmathoron's comments and also
|
Remove unnecessary argument.
|
@ingablundell Does your "thumbs up" still stand after the changes I made in response to @Silmathoron's comments? |
| @@ -337,7 +300,7 @@ class ht_neuron : public Archiving_Node | |||
| /** Timer (counter) for spike-activated repolarizing potassium current. | |||
| * Neuron is absolutely refractory during this period. | |||
| */ | |||
| int ref_steps_; | |||
| long ref_steps_; | |||
There was a problem hiding this comment.
I never thought about this before, but it could indeed be a good idea to generalize this change...
heplesser
changed the title
…tation in the second part.
| "- $m_{DK}^{\\infty}$ is a steep sigmoid which is almost 0 or 1 except for a narrow window around $d_{1/2}$.\n", | ||
| "- To the left of this window, $I_{DK}\\approx 0$.\n", | ||
| "- To the right of this window, $I_{DK}\\sim -(V-E_{DK})$.\n", | ||
| "The equations implemented in NEST are thus\n", |
There was a problem hiding this comment.
Either some other equation are missing or this is a wrong copy-paste
There was a problem hiding this comment.
Yes, that was a copy-paste leftover, fixed.
| " m_{DK}^{\\infty} &= \\frac{1}{1 + \\left(\\frac{d_{1/2}}{D}\\right)^{3.5}}\\\\\n", | ||
| " \\frac{dD}{dt} &= D_{\\text{influx}}(V) - \\frac{D-D_{\\text{eq}}}{\\tau_D} = \\frac{D_{\\infty}(V)-D}{\\tau_D} \\\\\n", | ||
| " D_{\\infty}(V) &= \\tau_D D_{\\text{influx}}(V) + {D_{\\text{eq}}}\\\\\n", | ||
| " D_{\\text{influx}} &= \\frac{D_{\\text{influx,peak}}}{1+ \\exp\\left(-\\frac{V-D_{\\theta}}{\\sigma_D}\\right)} \n", |
There was a problem hiding this comment.
Could you add in the "Note the following" list what you told me about the fact that, despite their ambiguous names, D_{influx} and D_{influx,peak}
represent ion concentration change
If you could elaborate a little on this here, I think it would be a good thing.
| "### Integration \n", | ||
| "\n", | ||
| "The original Synthesis implementation of the model uses Runge-Kutta integration with fixed 0.25 ms step size, and integrates channels dynamics first, followed by integration of membrane potential and threshold.\n", | ||
| "\n", |
There was a problem hiding this comment.
Maybe you could state which RK solver you are using?
| "\n", | ||
| "- The equation does not contain membrane capacitance. As a side-effect, all conductances are dimensionless.\n", | ||
| "- Na and K leak conductances $g_{\\text{NaL}}$ and $g_{\\text{KL}}$ are constant, although $g_{\\text{KL}}$ may be adjusted on slow time scales to mimic neuromodulatory effects.\n", | ||
| "- Reversal potentials are assumed constant.\n", |
There was a problem hiding this comment.
I would prefer: Reversal potentials
| "Synthesis: `ItChannel`\n", | ||
| "\n", | ||
| "##### Equations given in paper \n", | ||
| "\n", |
There was a problem hiding this comment.
Why are you calling it "equations" and not "definitions"?
There was a problem hiding this comment.
Because physicist are a sloppy bunch and call most things with a "=" in the middle an equation. Fixing this here may introduce inconsistencies elsewhere in the notebook, so I prefer to leave it as it is.
| "Intrinsic currents are based on the Hodgkin-Huxley description, i.e.,\n", | ||
| "\n", | ||
| "\\begin{align}\n", | ||
| "I_X &= g_{\\text{peak}, X} m_X(V, t)^N_X h_X(V, t)(V-E_X) \\\\\n", |
There was a problem hiding this comment.
I would prefer it if the meaning of the variables were stated explicitly. But maybe this is not necessary.
| "\n", | ||
| "#### Synaptic \"minis\"\n", | ||
| "\n", | ||
| "Synaptic \"minis\" due to spontaneous release of neurotransmitter quanta [HT05, p 1679] are not included in the NEST implementation of the Hill-Tononi model. This due to the fact that the total mini input rate for a cell was just 2 Hz and they cause PSP changes by $0.5 \\pm 0.25$mV only and thus should have minimal effect.\n", |
There was a problem hiding this comment.
Change that to a much simpler "..., because the total mini input rate".
| " t_{\\text{peak}} &= \\frac{\\tau_2 \\tau_1}{\\tau_2 - \\tau_1} \\ln\\frac{ \\tau_2}{\\tau_1}\n", | ||
| "\\end{align} \n", | ||
| "\n", | ||
| "where $t_{\\text{peak}}$ is the time of the conductance maximum and $\\tau_1$ and $\\tau_2$ are synaptic rise- and decay-time, respectively; $\\Theta(t)$ is the Heaviside step function. The equation is integrated using exact integration in Synthesis; in NEST, it is included in the ODE-system integrated using RK.\n", |
|
@ingablundell Thank you! I have revised the notebook according to your suggestions, with one exception. @apeyser Could you merge, as this is my PR? |
Implement the full NMDA model from the Hill-Tononi paper with instantaneous blocking and fast and slow unblocking. This PR will change model behavior.
As reviewers, I suggest @ingablundell and @suku248.